function predictions = predict_rsl(posterior, testdata);
%
% Performs predictions for test points based on the relevant subtask
% learning (RSL) model.
%
% Inputs: posterior = posterior approximation, obtained either from
% the initialization by initialize_rsl or after iterative optimization
% by optimize_rsl. testdata = matrix of test data; the model assumes
% that the test data come from the task of interest. 
% 
% The test data must be a similar matrix as the training data from the
% task of interest, so that each row is a data point, and in each row
% the first n_features elements are data features and the value of element 
% "n_features+1" is a constant one (1). There can be further elements after
% that but the prediction will not use them.
%
% Outputs: returns a vector of probabilities that the test points belong 
% to the second class (class +1), that is, p(class = +1|x) for each test
% point x.
% 


Sigma_U = posterior{1};
theta_U = posterior{2};

n_dim = length(theta_U);

X = testdata(:,1:n_dim)';
size(X);

logpredictions1 = (1 + 0.125*pi*sum(X.*(Sigma_U*X),1)).^0.5;
logpredictions2 = (theta_U'*X)./logpredictions1;

logpredictions = zeros(size(testdata,1),2);

% Predictions for class 1
logpredictions(:,1) = 1./(1+exp(-logpredictions2));
I = find(logpredictions(:,1) == 0);
I2 = find(logpredictions(:,1) > 0);
logpredictions(I2,1) = log(logpredictions(I2,1));
logpredictions(I,1) = logpredictions2(I);  % Approximation only!

if 0,
% Predictions for class -1
logpredictions(:,2) = 1./(1+exp(logpredictions2));
I = find(logpredictions(:,2) == 0);
I2 = find(logpredictions(:,2) > 0);
logpredictions(I2,2) = log(logpredictions(I2,2));
logpredictions(I,2) = logpredictions2(I);  % Approximation only!
end;

predictions = exp(logpredictions(:,1));

